//tihon.c

//regularization filter complex function magnitude computation
//by Tihonov method

#include <math.h>
#include "menu.h"
#include "user.h"
#include "message.h"
#include "chain.h"

static double flag=0,par,con,wl1,f01,n1,n2;
double besj(double);

//_________________________________________________________________________

void tihon_init(double f0, double wl, double a, double c,
		int m1,int m2,double size1,double size2,double step,
		double *min,double *max){

//********************** Parameter initialization **************************

double tihon(double,double);
double res;
float min1,max1;
       wl1=wl;
       f01=f0;
       n1=m1;
       n2=m2;
       con=c;
       par=(M_PI*a)/(wl*f0);
       flag=1;
       min1=tihon(0,0);
       max1=tihon(0,0);
//*    double u,v;
    double v;
    res=size1;
    if(size2>size1) res=size2;
 //   for(u=0;u<size1/2;u+=step){
       for(v=0;v<res/2;v+=step){
			    res=tihon(v,v);
			    if(res<min1) min1=res;
			    if(res>max1) max1=res;
    }
 /*   for(u=0,v=0;v<size2/2;v+=step*4){
			    res=tihon(u,v);
			    if(res<min1) min1=res;
			    if(res>max1) max1=res;
			    res=tihon(v,v);
			    if(res<min1) min1=res;
			    if(res>max1) max1=res;
    } */
 res=fabs(max1-min1);
*min=min1-0.05*res;

*max=max1+0.05*res;
return;
}
//___________________________________________________________________________

double tihon(double u,double v) {

//*********** Regularization filter magnitude at single point ***************

//*double h,h2,p,r,s,be;
double h,p,r,s,be;
      if(u==0 && v==0) {
	 u=0.0001;
	 v=0.0001;
      }
      r=sqrt(u*u+v*v);
      s=par*r;
      be=besj(s);
      h=2*be/s;
//      h2=fabs(h);
      p=h/( h*h + con * pow(u/(wl1*f01),n1) * pow(v/(wl1*f01),n2) );
      return  p; }


void filter_tikhonov(){
   char *name={"   DOE : Tikhonov filter   "};
   char *Par[]={
		  "Height             :",
		  "Width              :",
		  "Focal  Length      :",   //
		  "WaveLength (micron):",   //
		  "Discretization Step:",
		  "Defocus.  Parameter:",
		  "Noise/Signal Ratio :",
		  "First  Index       :",
		  "Second Index       :"   //
		 };
   char *Def[]={"    25.6          ",
		"    25.6          ",
		"    500          ",
		"    0.632         ",
		"    0.1           ",
		"     .2           ",
		"    0.0000001     ",
		"    3             ",
		"    3             "};

   sedit(25,5,9,Par,TmpPar,Def,name) ;

   if (!TmpPar[0]) return ;

   double  f, wl,  a, c;                 //
   int n1,n2;                                  //
//   double a_max, a_min;                        //

   sscanf(TmpPar[0],"%lf \0",&DoeHight);  sprintf(Def[0],"%.3lf",DoeHight);
   sscanf(TmpPar[1],"%lf \0",&DoeWidth);  sprintf(Def[1],"%.3lf",DoeWidth);
   sscanf(TmpPar[2],"%lf \0",&f);         sprintf(Def[2],"%lf",f);
   sscanf(TmpPar[3],"%lf \0",&wl);        sprintf(Def[3],"%lf",wl);
   sscanf(TmpPar[4],"%lf \0",&DoeStep);   sprintf(Def[4],"%lf",DoeStep);
   sscanf(TmpPar[5],"%lf \0",&a);         sprintf(Def[5],"%lf",a);
   sscanf(TmpPar[6],"%lf \0",&c);         sprintf(Def[6],"%.10lf",c);
   sscanf(TmpPar[7],"%d \0",&n1);         sprintf(Def[7],"%d",n1);
   sscanf(TmpPar[8],"%d \0",&n2);         sprintf(Def[8],"%d",n2);
   wl=wl/1000;
   if(MessageSize(DoeWidth,DoeHight,DoeStep)) return;

   tihon_init(f, wl, a , c , n1 , n2,DoeHight,DoeWidth,DoeStep,
				    &MinValueDoe,&MaxValueDoe);    // //

   chFun=chain(tihon," ");
   DoeDoe=Amplitude;              // //
   FormDoe=dftNonSymmetric;
   MessageName("Tikhonov Filter");

 //MessagePar(4,TmpStr);
   sprintf(TmpStr,"%s%s",Par[0],Def[0]); MessagePar(1,TmpStr);
   sprintf(TmpStr,"%s%s",Par[1],Def[1]); MessagePar(2,TmpStr);
   sprintf(TmpStr,"%s%s",Par[2],Def[2]); MessagePar(3,TmpStr);
   sprintf(TmpStr,"%s%s",Par[3],Def[3]); MessagePar(4,TmpStr);
   sprintf(TmpStr,"%s%s",Par[4],Def[4]); MessagePar(5,TmpStr);
   sprintf(TmpStr,"%s%s",Par[5],Def[5]); MessagePar(6,TmpStr);
   sprintf(TmpStr,"%s%s",Par[6],Def[6]); MessagePar(7,TmpStr);
   sprintf(TmpStr," N1=%d  N2=%d",n1,n2); MessagePar(8,TmpStr);
   MessageOut();

   UndoMenu(2) ;
   return;
}
